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I. INTRODUCTION 



Sulfidization impedes the use of palladium membranes for hydrogen sequester and cap- 
ture. Due to the presence of hydrogen sulfide (H 2 S) in any practical feed gas, sulfur rapidly 
covers the surface of the membrane rendering its catalytic properties inert. One possible 
way to overcome this difficulty is to alloy palladium with copper-^. In order to realize 
hydrogen capture and sequester, new membranes with different chemical compositions that 
resist sulfidization need to be devised, and their reactivity to sulfur needs to be examined. 
Computational methods offer an attractive way to rapidly prototype different chemical fam- 
ilies. 

In this paper we computationally model the temperature-dependent phase diagram for 
the Cu-Pd-S system using first principles methods. Prior experimental work- determined 
the ternary phase diagram for the Cu-Pd-S system. We seek a simple and semi-quantitative 
model to gain insight into the essential features of the system, in contrast to complicated 
models with many fitting parameters that can yield close agreement with experiment at the 
cost of obscuring the underlying physics. 

First, we describe our method for calculating ternary phase diagrams and apply it to Cu- 
Pd-S. Our model begins with first principles total energies, moves to focus on substitutional 
entropy, then examines vibrations. Comparison with experimental results shows general 
agreement. Certain details of energetics are explored through examination of the electronic 
density of states. We then apply our model to predict sulfidization threshholds for Pd 
interacting with S dimers and activities in Cu-Pd binaries. 

II. METHODS AND CALCULATIONS 

In order to calculate a phase diagram using first principles methods, all plausible crystal 
structures must be considered, in principle. In the case of Cu-Pd-S, a wealth of experimental 
crystallographic data gives an excellent starting point, especially for the binary compounds. 
We generalize these binaries to ternary structures formed from substituions of the ternary 
specie into either lattice sites or interstitial sites. Once the T=0K total energies for all likely 
crystalline structures have been calculated, it is straightforward to generate a T=0K phase 
diagram by creating the convex hull of total energies. 
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Subramanian and Laughlin^ determined the phase diagram for Cu-Pd experimentally, 
and Li et al. created a CALPHAD model 5 . Above T=871K, the entire composition range 
of solid Cu-Pd exists in fee solid solution. Below this temperature other solid phases are 
present. /3-CuPd has a CsCl structure and exhibits a Cu-rich composition range extending 
to a composition of 58% Cu. At copper concentrations near 75%, phases based on ID and 2D 
long period superstructures (LPS) compete. Studies using ab initio methods to examine this 
region of the phase diagram have already been performed^. As we are primarily concerned 
with the palladium-rich side of the phase diagram, we only calculate the ideal ID LPS, 
Cu3Pd.tP28, (using a notation of [chemical formula]. Pearson symbol) and ignore the 2D 
LPS and phase solubility ranges in the ID LPS. 

Determination of the phase diagram for Cu-S^ using first principles calculations is prob- 
lematic, as it contains an abundance of phases, many of which are non-stoichiometric or 
have large unit cells. We include only Cu-S structures with well-defined crystallographic 
refinements, as we do not expect this discrepency to affect our results on the palladium-rich 
side of the phase diagram. The Pd-S phase diagram^ is relatively simple, containing Pdi6S7, 
Pd4S, PdS, and PdS2 phases as line compounds at low temperature, with PdsS stabilized 
at T=829K. We summarize the phase diagrams here only for comparison; our method does 
not require pre-existing knowledge of the phase diagrams. 

In order to add temperature dependence to the phase diagrams, we introduce free energy 
models for select phases of interest, all of which are based on substitutional solid solution 
models with sublattice filling (where for some structures, the "sublattice" is the entire lat- 
tice). When calculating the ternary phase diagram and the activities, we work in the Gibbs 
ensemble, where we specify temperature, pressure, and chemical concentration. We neglect 
pressure dependence in the solid phases' free energies, as solid compressibilities are low. 
Phonon free energies are included, in the harmonic approximation, giving free energies of 
form 

Gtot G con fig Gphonon- (1) 

Based on our results, we chose to model only substitutional defects and not interstitials. 
We here summarize the configurational free energy G con fig of a simple substitutional solid 
solution model, valid in the dilute limit. In a substitutional solid solution, we have a set 
of site classes i, each with N; sites, an associated enthalpic cost of performing a chemical 
substitution of specie a into site class i, AE i Q , and a substitution concentration of x i Q ,. In 
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the substitutional solid solution model, it is assumed that each substitution is uncorrelated 
with all other chemical substitutions, and all spatial configurations have identical energies. 
This implies a linear dependence of AH with respect to composition, giving an enthalpic 
contribution of 

AH = AH + J2 x ha NAE ha) (2) 

i,a 

where AH Q denotes the enthalpy of the base structure with no substitutions. Including an 
ideal entropy of mixing 

AS = -k B ^ Ni[xi !a \og(xi, a ) + (1 - x it a) log(l - x i:(X )], (3) 

leads to a Gibbs free energy AG = AH — T * AS. Here, AHo and AE^ a can be determined 
from first-principles calculations. However, we control only the total impurity concentration 
for a species a, x Q , not the number of substitutions in a given site class x iQ . Since we 
work with a fixed number of sites, we fix the total concentration of substitutions, and thus 
minimize the free energy over all x i Q , subject to the constraint of specie number conservation. 

To better elucidate the model, we restrict to the special case of where only one type of 
substitution (with an enthalpic cost of AE) is allowed in only one sublattice containing N 
sites. A case with two sublattices will be considered later. The configurational free energy 
takes the form 

AG(N',T) = AH + AE * N' + kTN (x' hgx' + (1 - x') log(l - a;')) (4) 

where N' is the number of substitutions performed and x' = N'/N is the concentration of 
substitutions in the sublattice. The appropriate quantity to consider for thermodynamics is 
the free energy per atom (that is, the total number of atoms N), Ag = AG/N, expressed in 
terms of the x = N'/N, which can be shown to be 

Ag(x,T) = Ah + AE*x-TAS (5) 

where A/i = AH /N, and 7 = N /N is a measure of the number density of sublattice sites. 
Of particular importance is the entropic contribution to the free energy, 

— TAS = ksT'y^x log x + (7 — x) log(7 — x) — 7 log 7) (6) 

which we use for all phases modeled throughout this paper. 
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In order for a given substitution type to appreciably affect the free energy, the two 
necessary features are large 7 values, corresponding to a sublattice dense in the total lattice 
and thus a significant entropic contribution from the sublattice, and small positive (or any 
negative) enthalpic cost Ah for the substitutions, such that at a given temperature the 
entropic portion of the free energy can overcome the enthalpic cost (or reinforce the enthalpic 
gain) . 

We use VASP (the Vienna Ab-initio Simulation Package) — ^ a plane wave ab-initio 
package implementing PAW potentials^ to determine total energies. Previous work by 
Hu et al.— compared the accuracy of different pseudopotentials in the Pd-S system, and 
recommended PBEsol^ as most suitable for calculations involving sulfur. The total energy 
calculations were performed by fully relaxing atomic positions and lattice parameters until 
energies converged to within 0.1 meV/atom. Subsequent convergence in k-point density was 
performed. A common energy cutoff of 273 eV was used. In order to model solubility ranges, 
calculations were performed both in unit cells and larger supercells. All phonon calculations 
were done using density functional perturbation theory, and were performed in at least a 
2x2x2 supercell of the unit cell, with the exception of PdigS7.cI46, where due to its large 
size only a 2x2x2 supercell of the primitive cell was used. Table [J shows a list of structures 
and their total energies. 



III. RESULTS AND DISCUSSION 
A. Phase Diagrams 

Given total energies, we calculate the enthalpy of formation for our structures at T=0K, 



where Ah is the calculated T=0K total energy for the base structure, Xj is the composition 
of a given specie i in the structure, and Ahi is the T=0K total energy for pure specie i 
in a chosen reference structure. By convention, this reference structure is chosen to be the 
equilibrium structure at T=0K, so that at T=0K all enthalpies of formation for the pure 




(7) 
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Phase 


Pearson Symbol 


Prototype 


Group # 


Space Group 


AH 


dE 




cF4 


Cu 


225 


Fm3m 








pd 16 


cF4 


Cu 


225 


Fm3m 








S2 Dimer 


N/A 


N/A 


N/A 


N/A 








Cu 3 PdiZ 


tP28 


Cu 3 Pd 


99 


P4mm 


-103.0 





CuPd 


cP2 


CICs 


221 


Pm3m 


-116.9 





Cu 2 S^ 


mP144 


Cu2S 


14 


P2i/c 


-496.4 





CuS^ 


oC24 


CuS 


63 


Cmcm 


-715.4 





CuS^ 


cP12 


FeS 2 


205 


Pa3 


-821.5 





Pd 4 S^ 


tPIO 


Pd 4 Se 


114 


P42iC 


-418.9 





Pd^S^ 


cI46 


Pdi 6 S 7 


217 


I43m 


-589.7 





PdS^i 


tP16 


PdS 


84 


P4 2 /m 


-892.8 





PdS# 


oP12 


PdSe 2 


61 


Pbca 


-952.4 







hP12 


CuS 


120 


I4c2 


-712.9 


2.5 


Cu 2 S^ 


tP12 


Cu 2 S 


96 


P4 3 2 4 2 


-479.3 


3.3 


Cu 4 Pd21 


cP4 


AuCu 3 


221 


Pm3m 


-99.7 


3.4 


Pd 3 S^ 


0CI6 


Pd 3 S 


40 


Ama2 


-492.8 


8 


CusoPdscri 


cF4 


Cu 


225 


Fm3m 


-106.9 


10 


Cu 4 Pd^ 


tP20 


Cu4Pd 


84 


P4 2 /m 


-69.6 


12.8 


Cu 2 Pd 3 S 4 


mP18 


Cu 2 Pd 3 Se 4 


14 


P2i/c 


-743.0 


13.2 




cF12 


CaF 2 


225 


Fm3m 


-503.1 


20.6 


Cu 


cI2 


W 


229 


Im3m 


36.3 


36.3 



TABLE I: A partial list of calculated structures, focusing on "base" structures with no 
substitutions. All energetic quantities are given in units of meV/atom. Cu 2 Pd 3 S 4 is a 
hypothetical structure based on Cu 2 Pd 3 Se 4 .mP1822,. Cu.cI2 and CuPd.cP2 are 
hypothetical structures used to model the /3-CuPd phase. AH and dE are defined in eqs. [7] 

and[8] 
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species vanish. In the Cu-Pd-S system, the standard reference structures would be Cu.cF4 
for copper, Pd.cF4 for palladium, and S.mP48 for sulfur. 

As we are interested in interaction of gaseous sulfur with copper-palladium membranes, 
we have chosen to use gaseous sulfur dimers rather than the standard S.mP48 when calcu- 
lating phase diagrams. In the case of the ternary phase diagram, all structures of interest 
have relatively low sulfur composition, so the mismatch in reference phase compared to 
experimental phase diagrams is negligible. 

The quantity of importance for constructing temperature dependent phase diagrams is 

dE = Ah-[Ah], (8) 

where [Ah] is the enthalpy of formation of the convex hull at a structure's stoichiometry. 
By definition, this quantity is non-negative, and it vanishes only when the structure lies on 
the convex hull. This quantity is a measure of how far above the convex hull a particular 
structure lies at T=0K. Structures with large dE's require large entropic contributions to 
free energy in order to overcome enthalpic costs. 

Figure [1] shows the T=0K binary enthalpies of formation. Observe that palladium and 
sulfur react strongly, as can be seen by the maximum enthalpy of formation of 0.95 eV/atom. 
This strong enthalpy of formation reflects the tendency of pure palladium membranes to 
sulfidize. Next, copper and sulfur react almost as strongly, with a maximum enthalpy of 
formation of 0.82 eV/atom. Finally, copper and palladium have a maximum enthalpy of 
formation of only 0.12 eV/atom, indictating that the interspecies binding between copper 
and palladium is relatively weak. This is expected from the Cu-Pd phase diagram, as all 
known phases are variants of fee (the T=0, p=0 structure for both Cu and Pd) or bec 
(which is reached from fee by a martensic transition^). We find that the known structure 
CuS.oC24 is mechanically unstable, with an imaginary phonon mode that stabilizes a slight 
monoclinic distortion, and this distorted monoclinic structure, CuS.mC24, is the true ground 
state structure according to DFT. The stable phases predicted correspond to known low 
temperature stable phases, with the exception of CuS2-tP12, which is not observed stable at 
temperatures as low as 300K (its stability is mostly likely due to the use of S2 as a reference 
structure), and CuS.mC24, the experimentally observed stable structure being CuS.hP12. 

Experimental results^ show Pdi6S7.cI46 and Pd4S.tP10 having solubility ranges for cop- 
per, which implies the existence of substitutional or interstitial defects. Pdi 6 S 7 .cI46 contains 
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FIG. 1: (color online) Binary enthalpies of formation for Pd-S, Cu-S, and Cu-Pd. The x 
axis corresponds to concentration, and y axis corresponds to enthalpy of formation. Black 
symbols denote stable structures, blue correspond to structures witin 3 meV of the convex 

hull, and red to structures above 3 meV. Circles with thick borders are known low 
temperature stable structures, circles with thin borders are known high temperature stable 
structures, and squares denotes all other cases (including hypothetical structures). For 
Cu-Pd, to better distinguish phases, filled symbols correpond to bcc structures, empty 
symbols to fee, and dashed for the three other structures (CusPd.tP28, Cu4Pd.tP20, and 

the SQS configuration for fee at 50% Pd). 






four Wyckoff site classes: 8c and 24g sites containing palladium and 8c and 6b sites con- 
taining sulfur (see Fig. [2]). For Pd4S.tP10 there are only two site classes, an 8g containing 
palladium and 2a containing sulfur. 

Shown in Figure [3] is our calculated T=0K phase diagram for the Cu-Pd-S system (sup- 
porting data is in Table [TT]) . We find the 8c site class in Pd^Sr is favorable to copper 
substitution. At T=0K it is enthalpically favorable to occupy 6 of the 8 possible lattice 
sites, with a sudden upturn for further filling of the 8c site class. The only other substi- 
tutions that are found favorable are S and Cu in /3-CuPd, Cu in Pd.cF4, and Pd in CuS 2 . 
Substitutions of Cu in the /3-CuPd phase have a solubility range from 50% Cu to 55.5% at 
T=0K. We find the /3-CuPd phase is only favorable to Cu substitutions on the Pd sublat- 
tice, and not Pd substitutions on the Cu sublattice, which is experimentally supported by 
observations that the palladium content of /3-CuPd never exceeds 50%. As none of the Cu- 
and Pd- rich structures are in competition with CuS2-tP12, and the entropic contribution is 
small, we ignore substitutions in CuS2-tP12. We also find a small sulfur solubility range in 
the Cu sublattice of /3-CuPd at T=0K, which is supported by experimental phase diagrams 
at T=673K and T=823K. 

While these are the only stable T=0K substitutions, we expect that at higher tempera- 
tures, substitutions into other sites that had suitably low enthalpic cost will be stabilized 
by entropy: PdS.tP16 (2c), Pdi 6 S 7 .cI46 (24g), and Pd 4 S.tP10 (8e) sites being possible can- 
didates (see Table ITT]) . However, we chose to exclude modeling the PdS.tP16 2c substitution 
range, as due to the small concentration of 2c sites in PdS.tP16, combined with competi- 
tion with other phases with strong entropic effects, we expect the concentration range to be 
negligible. 

To determine solubility ranges more precisely, we model phases analytically (see Table 
IHIj) . using first principles calculation data to suggest suitable enthalpic models. We choose 
the models to smooth out small flucuations in the first principles calculated energies while 
still keeping the essential trends seen intact. The phases modeled are Pd4S, Pd^Sy, /3-CuPd, 
and the CuPd-FCC phase. In all cases we use the ideal entropy (Eq. E]) for sublattice filling. 
For Pdi 6 S 7 , we use a piecewise linear function (see Fig. [H]) for the enthalpic part of the free 
energy, because the enthalpic contribution to the free energy is increasingly favorable up to 
6 substitutions (x=0.13) on the 8c sublattice, but from 6 to 8 substitutions (x=0.174) the 
favorability is reduced. Previous crystallographic work by Matkovic 23 observed a crystalline 
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FIG. 2: (color online) Crystal structures for P 4 S.tP10 and P16S7.CI46. The unit cells are 
outlined. Atomic sizes indicate vertical height. 
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phase of P0I16S7 with 14% Cu concentration, the Cu occupying 8c sites. However this work 
was performed at 853K and did not suggest any mechanism by which these substitutions 
are stabilized. For Pd4S, it was necessary to add a quadratic term to the free energy to 
accurately fit our first principle results. Both CuPd-FCC and /3-CuPd were fit using a 
first-order regular solution model, with CuPd-FCC having the entire lattice available for 
substitution and /3-CuPd having only the Pd sublattice available for substitution by Cu. 
We also compute a SQS configuration for fee at 50 Pd% 34 with purely random correlations 
until the 8th nearest neighbor and find it lies above the convex hull by 40 meV/atom. 

Experimental phase diagrams show a sulfur solubility range near Cuo.sPdo.s trending 
in the palladium-rich direction. We find dilute S substitutions onto the Cu sublattice in 
/3-CuPd, creating a phase trending in the palladium-rich direction, are stable at T=0K. 
Subsequent substitutions become rapidly unfavorable. From this, we expect limited entropic 
stabilization of additional substitutions of S in /3-CuPd at higher temperatures, leading to 
a roughly constant solubility range. 

To include phonon free energy in solid solution phases, we calculate phonon free energies 
in the harmonic approximation at well-ordered configurations at select compositions and 
linearly interpolate between these compositions across the phase region. To model the free 
energy of gaseous sulfur, we use emperical data taken from the NIST-JANAF tables, as 
explained later in the paper. 



Figures H] and [5] compare our calculated phase diagrams with experiment for T=673K and 
823K, respectively. Our T=673K phase diagram matches the experimental results qualita- 
tively. Note, however, that the "Cu2_ x S" region of the experimental phase diagram contains 
multiple structures, many of which are non-stoichoimetric, whereas ours contains only one 
structure. The phase labeled "CusPd" (or ID LPS) has a slight solubility range experimen- 
tally whereas we model it as a well-ordered line compound. Our T=823K diagrams show 
more disagreement, as there is only one Cu-S structure that is stable and a PdsS structure 
is stabilized at this temperature. More recent phase diagrams have PdsS being stabilized 
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FIG. 3: (color online) The calculated T=0K phase diagram for the Cu-Pd-S system. 

Plotting symbols as in Fig. [TJ 



at T=829K, however. This is likely due to entropic effects other than configurational and 
phononic (in the harmonic approximation) that were not included. As can be seen, PdaS 
is part of our collection of calculated structures but is not a stable low-temperature phase. 
At both temperatures we attempt to model sulfur solubility in the /3-CuPd phase but find 
sulfur solubility to be unstable at temperatures of interest. This is surprising, as our T=0K 
phase diagram has a sulfur solubility range in /3-CuPd similar in shape to what is experimen- 
tally observed at higher temperatures, and the enthalpic behavior of our sulfur substitutions 
supports the observed temperature independence of the sulfur solubility range. 

The Pdi 6 S7.cI46 copper solubility range is not appreciably affected by the temperature, 
as it is dominated by enthalpy rather than entropy. The 24g site is unfavorable for copper 
substitution. The high multiplicity of the 24g site class gives a large entropic contribution 
to the free energy at high temperature relative to the 8c site class. However, we find that 
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Base Structure 


Site Class 


Chemical Specie 
Change 


Enthalpic Cost 
(meV/atom) 


Maximum Impurity Concentration 
in Site Class 


/3-CuPd 


2a 


Pd to Cu 





12.5 


/3-CuPd 


2a 


Cu to S 





3.8 


/3-CuPd 


2a 


Pd to S 





3.8 


Pd.cF4 


4a 


Cu to Pd 





21.9 


Pdi 6 S 7 .cI46 


8c 


Pd to Cu 





12.5 


Pdi 6 S 7 .cI46 


24g 


Pd to Cu 


12.9 


8.33 


Pd 4 S.tP10 


8e 


Pd to Cu 


4.7 


3.13 


CuS 2 .cP12 


4a 


Cu to Pd 





6.25 


PdS.tP16 


2c 


Pd to Cu 


5.3 


50.00 



TABLE II: A list of substitutions into binary structures to produce ternary structures. 
Only low enthalpic cost substitutions are shown. Impurity concentration denotes the 
concentration in the appropriate sublattice, not concentration in the structure as a whole. 



the enthalpic cost (see Figure [6]) associated with occupying 24g sites is too high relative to 
the entropic contributions at the given temperatures and only affects the solubility range 
by 1% to 2%. Thus the solubility range consists almost completely of 8c sites and is unap- 
preciably affected by temperature. For this reason, our model (Table IIIII) only includes 8c 
sites and not 24g sites. For Pdi6S7.cI46, at both measured temperatures the experimental 
copper solubility limit is 20%, and our calculated value is 17.5%, the mismatch likely due 
to Cu occupation of 24g sites. Our results explain why this solubility range appears to be 
independent of temperature; it is not stabilized due to entropic contributions but due to 
enthalpy arising from details of the electronic structure (see subsection IIIIBI) . 

In contrast, Cu substitutions in Pd4S.tP10 are not enthalpically favorable at T=0K but 
are entropically stabilized at higher temperature. As can be seen in the above phase dia- 
grams, the copper solubility for Pd4S.tP10 is temperature-dependent in both the calculated 
and experimental work. While our results agree qualitatively with experimental results, 
there is still a discrepancy in the copper solubility range. For Pd4S.tP10, at T=673K the 
experimental value is 12.7% and ours is 2.5%, and for T=823K the values are 18.6% and 5% 
respectively. Solid solution models overestimate the configurational entropy, as including 
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Phase 


Substitution Type 
(Cu Composition Range) 


7 


Phonon Concentrations 


hconfig 

(eV/atom) 


Pd 4 S 


Cu@Pd8e (0-0.8) 


8/10 


0, 0.8, 


-5.350 + 1.612x + 0.225x 2 


Pdi 6 S 7 


Cu@Pd8c (0-0.175) 


8/46 


0, 0.175 


-5.237 + 1.445x(x < 6/46), -5.251 + 1.554x(x >= 6/46) 


/3-CuPd 


Cu@Pd2a (0.5 - 1) 


1/2 


0.5, 1 


-5.744 + 1.773a; + (x - 0.5)(1 - x)(-0.206(x - 0.5) - 1.526(1 - x)) 


CuPd-FCC 


Cu@Pd4a (0-1) 


1 


0, 0.50, 0.75, 1 


-5.477 + 1.469x + x(l - x)(-0.529x - 0.267(1 - x)) 



TABLE III: Description of our phase models. 7 denotes the number concentration of the sublattice in the total lattice and is 
used for scaling the entropic contribution to energy (see equation |6]). "Phonon concentrations" denotes the concentrations for 

which phonon free energies were calculated, x denotes the copper composition. 



QU Cu 3 Pd B-CuPd p^J 

(a) Experimental, T=673K 




Cu 0.9 0.8 0.7 0.6 0.5 0.6 0.7 0.8 0.9 Pd 

Cu 3 Pd CuPd-FCC 



(b) Calculated, T=673K 

FIG. 4: (color online) The experimental phase diagram for the Cu-Pd-S system at 
T=673K, and our calculated phase diagram. 

correlations between substitutions will alter the energies of distinct configurations with the 
same composition and reduce the multiplicity. We have thus overestimated the entropy of 
the Pdi 6 S 7 phase, whose solubility range extends to almost complete filling of the 8c site, 
relative to the Pd 4 S phase, whose filling in the 8e sublattice is sparse. Thus we would expect 
more accurate free energy models will lead to decreased entropy in the PdigSy phase with 
comparatively small changes in Pd4S, increasing the solubility range of Pd4S. This will not 
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(b) Calculated, T=823K 

FIG. 5: (color online) The experimental phase diagram for the Cu-Pd-S system at 
T=823K, and our calculated phase diagram. 



affect the stability of PdigSy as it is primarily enthalpically stabilized. 

The large solubility range in the palladium rich side of the CuPd-FCC phase at T=0K is 
also a factor in the incorrect solubility ranges for Pd4S.tP10, as we have found that altering 
the fitting parameters for the CuPd-FCC phase appreciably affects the solubility range for 
Pd 4 S.tP10. As we also encountered an issue with sulfur solubility in /3-CuPd, with the sulfur 
solubility range destabilizing relative to CuPd-FCC at temperature below temperatures of 

16 




FIG. 6: (color online) Comparison for the enthalpy of formations for filling 8c sites versus 

filling 24g sites in Pdi 6 S 7 .cI46. 



interest, and we have a large solubility range on the Pd-rich side at T=0K for the CuPd-FCC 
phase, it appears that one major source of error is inaccurate modeling of the CuPd-FCC 
phase. We will return to this issue in the activities portion of this paper. 

For our Cu-Pd binary at high temperatures, we find that Cu3Pd.tP28 becomes unstable 
relative to the CuPd-FCC phase at T=827K, and /3-CuPd becomes unstable relative to 
CuPd-FCC and Cu3Pd.tP28 at T=536K and a composition of 59.5%, significantly below the 
experimentally observed temperature of T=871K but close to the experimental composition 
of 58%. This relative instability of /3-CuPd at temperatures of interest is responsible for the 
observed instability of any sulfur solubility in /3-CuPd in our T=673K and T=823K phase 
diagrams. 

But even with these discrepancies, it is clear that we have predicted major details of the 
Cu-Pd-S ternary phase diagram, specifically the mechanism for stabilization of the solubility 
ranges of Pd 4 S.tP10 and Pdi 6 S 7 as well as details about the stability of the /3-CuPd system, 
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E [eV] 

FIG. 7: (color online) Calculated electronic densities of states for Pdi 6 S7.cI46 with copper 
substitutions at the palladium 8c sites. Dashed lines indicate the Fermi level Ep. The inset 
shows the DOS relative to Ep and enlarges the pseudogap region. 

using only first principles calculations with no empirical fitting. 

B. Electronic Structure 

The question remains why the Pdi6S7.cI46 8c site class favors copper in the first 6 substi- 
tutions, but substitution in the final two 8c sites are unfavorable. Consider the constituent 
binaries: copper binds less strongly with sulfur than palladium, and copper binding with 
palladium is relatively weak, so we would expect copper substitution into Pd-S binaries to 
be slightly unfavorable. If it were to turn out to be favored, then we would expect Cu to 
fully occupy the 8c class. First-principles methods provide electronic structure information 
that can resolve this issue. 

Figure [7] shows the electronic densities of states for substitutions in Pdi 6 S7.cI46. Each 
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colored solid line corresponds to the electronic density of states for a different number of 
substitutions in the 8c site class of Pdi6S7.cI46, and the colored dashed lines correspond to 
the associated Fermi levels. It can be seen that the densities of states follow a rigid band 
model. Here copper is being substituted for palladium, and copper has one more electron in 
its valence shell than palladium. Under a rigid band model, performing a substitution adds 
one electron per unit cell without appreciably changing the valence electronic eigenstates, 
increasing the Fermi level relative to the fixed band structure. 

In all pseudogap lies near the Fermi level. In the case of no substitutions, 

the pseudogap lies to the right of the Fermi level. As the number of copper substitutions 
increases, the increased electron count drive the pseudogap to the left, towards the Fermi 
level. The inset shows the densities of states relative to the Fermi levels. At six substitutions 
the Fermi level lies almost directly at the minimum of the pseudogap. This is a well-known 
stabilization mechanism for alloys^. As the number of substitutions increases above six, the 
Fermi level is driven away from the pseudogap, destabilizing the structure and leading to 
the sudden increase in the enthalpy of formation (see Fig. |6]). This pseudogap stabilization 
mechanism at T=0K, combined with large enthalpic cost of 24g substitutions, explains the 
temperature independence of the Pdi 6 S 7 solubility range. 

C. Predominance Diagrams 

We now calculate the predominance diagram for the Pd-S2 system. We work with a 
fixed quantity of palladium interacting with a reservoir of gaseous S2, with temperature and 
pressure being free parameters. For an ideal gas, this implies that we control the chemical 
potential of the gas, and thus we work in the semi-grand canonical ensemble: 

S(T,p,x) =G(T,p,x) - /j,s 2 N S2 (9) 

where G is the previously defined Gibbs free energy, /ig 2 is the chemical potential for the 
sulfur dimers, and Ng 2 = N$/2 is the half the number of sulfur atoms in a given phase 
(as the chemical potential is for sulfur dimers). To determine sulfidization threshholds, we 
calculate the semi-grand canonical free energy for all phases of interest, Pd.cF4, Pd4S.tP10, 
Pdi6Sy.cI46, Pd3S.oC16, and PdS.tP16. The phase with the lowest free energy will be stable. 
All Pd-containing phases of interest are line compounds, containing only total energy and 
phonon contributions. 
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The chemical potential for S2 complicates the calculation, as it is a gaseous phase. We in- 
corporate empirical results into our method, in a manner consistent with the results obtained 
from first principles calculations. For the free energy of the gas, we have 

n(T,p) = n°(T) + kTlog?-. (10) 

Po 

Here \x is the chemical potential per molecule, po is a reference pressure (by convention taken 
to be 1 bar, as will be done here), and /jP is the chemical potential per molecule measured 
at the reference pressure, 

fi°(T) = h-Ts \ p=1 , (11) 

where H and S are obtained from the NIST-JANAF tables 36 . /i° obtained from the JANAF 
tables has a physically unrealistic infinite value at T=0K due to the presence of 1 /T terms in 
both the H and TS contributions. As shown in Figure [8], /x° is nearly linear for temperatures 
of interest, so we fit to a linear equation using the 770K < T < 1430i^ region to eliminate 
the T=0K pole. There is still an issue of constant offsets to energy; the JANAF tabulated 
values are relative to the enthalpy at T=298K, and whereas all our calculated values are 
computed at T=0K. In order to incorporate the empirical data into our method, we must 
shift its zero of energy. In particular, at T=0K, we set the empirical data to match what 
first principles calculations would have predicted, 

gS 2 \t=0K— & con fig + ^phonon \t=0K ■ (12) 

As we have measured all energies relative to the tie-plane between pure elements, e con fi g = 0, 
and so the zero-point energy e p h onon is the only contribution at T=0K taken from first 
principles calculations. 

Shown in Table IIVI are our predicted sulfidization threshholds relative to experimental 
predominance diagrams obtained by Taylor—. Predictions with and without phonon free 
energies are shown. It is necessary to include phonons to obtain reasonable agreement 
with experimental results, as when phonon free energies are excluded, the relative errors in 
pressure are as large as four orders of magnitude. When including phonon free energies, in 
all cases our pressure estimates are within a factor of 6 of experimental results. It should also 
be be noted that Taylor's phase boundaries are themselves based on fits from experimental 
enthalpy data. One difficulty is that Taylor has a liquid region in the high-temperature, low 
pressure region of the predominance diagram, whereas we are only interested in solid phases, 
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FIG. 8: (color online) The chemical potential of S2 dimers. Black is NIST- JANAF, and red 

is our linear fit. 



so experimental threshholds from the solid regions were extrapolated to liquid regions to give 
numbers suitable for comparison. This extrapolation is valid as Taylor's threshholds in the 
solid regions were found to be extremely linear. Another difficulty is that PdaS is stable for 
a small pressure range in between temperatures of 830K and 920K, and it would likely be 
stable for higher temperatures if not for the liquid phase, but our calculations do not predict 
it stable in this region. This is most likely due to errors in the phonon free energy, which as 
already shown appreciably affects results. 



D. Activity Calculations 

In this last section, the activities of species in binary Cu-Pd are calculated, as these can 
be used to predict sulfidization threshholds for Cu-Pd in the presence of S. We wish here to 

21 



T (K) 


Transition 


logiQ\fP Calculated 
Without Phonons 


logio\^P Calculated 
With Phonons 


logiQ^/P Experimental^! 


1000 


PdS to Pd 16 S 7 


-1.1 


-1.8 


-2.2* 


1000 


Pdi 6 S 7 to Pd 4 S 


-1.4 


-2.6 


-2.4* 


1000 


Pd 4 S to Pd 


-3.3 


-5.3 


-5.2 


800 


PdS to Pd 16 S 7 


-3.1 


-3.5 


-3.7 


800 


Pd 16 S 7 to Pd 4 S 


-3.5 


-4.5 


-4.7 


800 


Pd 4 S to Pd 


-5.9 


-7.6 


** 



TABLE IV: Sulfidization threshhold pressures for the Pd-S2 system. Calculated values 
with and without phonon contributions are seperately shown. Values * are liquid 
experimentally, but were extrapolated from the solid region of the predominance diagram 
for comparison purposes. ** not indicated in experimental results. 



find the activities across the entire composition range with multiple phases. 

The definition for the activity of a specie A as a function of composition x A , pressure P, 
and temperature T is^ 

H A {x A ,P,T)-fi° A {P,T) 



log a A (x A ,P,T) 



where a A is the activity for specie A, fi A (x A ) is the chemical potential for specie A, and fi A 
is a reference chemical potential for specie A. As we are specifying x A , P, and T, the natural 
choice of thermodynamic potential is the Gibbs free energy, 

G = G{N A ,N B ,p,T) = haNa + VbN b , 

where N A and N B are the total number of A and B atoms, respectively, and the thermody- 
namic derivative giving the chemical potential is 

dG 



Va(Na,N b ,P,T) 



ON, 



However, the quantity of importance for phase transitions is the Gibbs free energy per atom, 

G(N A ,N B ,P,T) 
g = g(x A ,P,T)= Na + Nb ■ 

We may express the chemical potential in terms of intensive quantities as 



Ha(x a , P,T)=g + (l- x A ) 



dg 

dx A 



\P,Ti 
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where X{ is the composition of species i, and we have chosen to parameterize the free energy 
by the composition for A. Similarly, 

Hb(x a ,P,T) = g + (1 -x B )-^- \ pt = g - x A -^- | PT . 

ax B oxa 

These activities depend only on the free energy and the derivative of the free energy at a given 
composition. They are piecewise analytic functions across the composition range, completely 
determined by the thermodynamically stable phase at that composition. The activities can 
therefore be viewed as alternative representation of the underlying phase diagram for the 
system. They are discontinuous only when the derivative of the free energy changes (the 
free energy itself must always be continuous, due to the requirement of convexity). 

As mentioned before, we believe the magnitude of the calculated free energy of our fee 
phase is artificially large, destabilizing the (3 phase, and thus any calculated activities would 
be flawed near the 50% Cu region of the phase diagram. To support this assertion, here 
we refit the fee phase using only low-solubility first principles data on either side of the 
composition range. This reduces the magnitude of the free energy of the fee phase while still 
using first principles derived calculations. This should not be viewed as abandoning our pure 
first principles phase diagram as shown before, but rather a correction inspired by known 
experimental results presented in an alternative-but-equivalent format. The SQS estimate 
for fee at 50% Pd lies almost directly on this free energy curve, further supporting its usage 
for quantitative estimates. We also neglect phonon free energy for activity calculations, as 
for ordered fee and bec at 50% Pd, the difference in phonon free energy per atom is only 2 
meV at T=1000K. 

Shown below in Figure [9] in solid lines are our calculated activities for palladium in the 
Cu-Pd binary, using the free energy models from Table the modified fee phase, and our 
calculated data for CusPd. The /3 phase is the sharply increasing region in the center of 
the graph. This sharp behavior is due to the entropic portion of the free energy, which has 
a logrithmic singularity at x=0.5 and sharply varying derivitive near this singularity. The 
smooth curves on the left and right sides of the graph for T<1000K denote the fee phase. 
Linear free energy functions give constant activities as a function of composition, and thus 
phase coexistence regimes appear as plateaus in activities, at a fixed temperature. Two sets 
of plateaus are present. The set of plateaus on the left side are coexistence regimes between 
the bec and fee phases. There are 3 types of plateaus on the right side. For T < 900K, the 

23 



plateaus are due to coexistence with the line compound CusPd, where the plateaus in the 
region 0.6 < x < 0.8 indicate coexistence with bcc and the plateaus at 0.8 < x < 0.9 indicate 
coexistence with fee. The plateau on the right side of the plot for T=900K is a special case. 
Our free energy models predict that the C^Pd phase becomes thermodynamically unstable 
at 899K and fee becomes stable, so the plateau at 900K indicates coexistence between bcc 
and fee. At a temperature of 960K, our model predicts that the bcc phase becomes unstable 
relative to the fee phase, disappearing at a composition of x=0.576. By decreasing the fee 
free energy magnitude we have properly stabilized the (3 phase at temperatures of interest 
and brought the transition temperature to within 90K of the expected value. The 1000K 
and 1350K plots are pure fee throughout the entire composition range. 

The activities of Pd at various temperatures were also calculated using CALPHAD 
(acronym of CALculation of PHAse Diagrams) method. The thermodynamic database 
developed by Li et al 5 was used. Our activities are systematically low compared to the 
CALPHAD activities (shown in dashed lines on Figure [9]), most likely due to our ideal en- 
tropy approximation. As the true entropy in solid solutions is reduced relative to the ideal 
entropy, the true free energy will be greater, leading to the true activities being greater than 
our calculated activities. This is especially true for fee, where the refitted free energy still 
overestimates the magnitude of the free energy and leads to qualitatively different behavior 
in the high Pd limit. CALPHAD results give activities that exceed Raoult's Law in the high 
Pd limit, which would require positive enthalpies of formation for Cu substitution in Pd 
using a solid substitution model in the dilute limit. This contradicts first principle calcula- 
tions which unambiguously predict Cu substitutions in bulk Pd having negative enthalpies 
of formation (see Figure [TJ. Clearly more precise models are needed to understand the 
binary in the high Pd limit. However, for the (3 phase and associated coexistence regimes, 
the first principles and CALPHAD results are in close agreement, and if one is interested 
in regions outside the single-phase fee region, first principles calculations gives reasonably 
accurate results even with simple models. 

IV. CONCLUSIONS 

Using only first principles calculations and solid solution models, we were able to semi- 
quant it at ively model the ternary phase diagram for the Cu-Pd-S system without any em- 
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FIG. 9: (color online) Activities of Pd in the Cu-Pd binary. First principles results are 
denoted with solid lines, and CALPHAD with dashed lines. 



pirical fitting parameters. We obtained the sulfidization threshholds for the Pd-S2 system 
and activities that reflect the essential features of the Cu-Pd binary system at tempera- 
tures where multiple phases exist across the composition range. More detailed models are 
necessary to reach better quantitative agreement with experimental results. However, our 
work reveals the essential physics underlying the phase behavior. For example, we distin- 
guish enthalpic vs. entropic driven substitution, and we identify a pseudogap stabilization 
mechanism, yielding a temperature-independent solubility range for Cu in Pdi 6 S7. 
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